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Abstract 

This short note presents the Lambert W(x) function and its possible application in 
the framework of physics related to the Pierre Auger Observatory. The actual numerical 
implementation in C++ consists of Halley's and Fritsch's iteration with branch-point 
expansion, asymptotic series and rational fits as initial approximations. 
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Figure 1: The two branches of the Lambert W function, W_i(jt:) in blue and Wo(x) 
in red. The branching point at (— e^^, —1) is denoted with a green dash. 

1 Introduction 

The Lambert W(x) function is defined as the inverse function of 

yex^py = x, (1) 

the solution being given by 

y = w(x), (2) 

or shortly 

W(x)expW(x) = X. (3) 

Since the x ^ x exp x mapping is not injective, no unique inverse of the x exp x function 
exists. As can be seen in Fig. 1, the Lambert function has two real branches with a branching 
point located at (— The bottom branch, W_i(x), is defined in the interval x G 
[— e^^, 0] and has a negative singularity for x — )■ 0^. The upper branch is defined for 
X G [—e^^, oo]. 

The earliest mention of problem of Eq. (1) is attributed to Euler. However, Euler himself 
credited Lambert for his previous work in this subject. The W(x) function started to be 
named after Lambert only recently, in the last 10 years or so. The letter W was chosen by 
the first implementation of the W(x) function in the Maple computer software. 

Recently, the W(x) function amassed quite a following in the mathematical community. 
Its most faithful proponents are suggesting to elevate it among the present set of elementary 
functions, such as sin(x), cos(x), ln(x), etc. The main argument for doing so is the fact that 
it is the root of the simplest exponential polynomial function. 

While the Lambert W function is simply called W in the mathematics software tool 
Maple, in the Mathematica computer algebra framework this function is implemented under 
the name ProductLog (in the recent versions an alias LambertW is also supported). 

There are numerous, well documented applications of W(x) in mathematics, physics, 
and computer science [1, 3]. Here we will give two examples that arise from the physics 
related to the Pierre Auger Observatory. 
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Figure 2: Left: The Moyal function M(z). Right: A family of one-parametric Gaisser- 
Hillas functions g{x; Xmax) for x^ax in the range from 1 to 10 with step 1. 



1.1 Moyal function 

Moyal function is defined as 

M(x) = exp (-^ (x + exp(-x))) . (4) 
Its inverse can be written in terms of the two branches of the Lambert W function, 

M±^(x) = Wo,_i(-x^) -21nx. (5) 

and can be seen in Fig. 2 (left). 

Within the event reconstruction of the data taken by the Pierre Auger Observatory, the 
Moyal function is used for phenomenological recovery of the saturated signals from the 
photomultipliers. 



1.2 Gaisser-Hillas function 

In astrophysics the Gaisser-Hillas function is used to model the longitudinal particle density 
in a cosmic-ray air showers [4]. We can show that the inverse of the three-parametric 
Gaisser-Hillas function. 
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is intimately related to the Lambert W function. Using rescale substitutions, 

X — Xo , 

X = 7 and 

A 

_ Xmax ~ Xq 
^max — ~^ / 

the Gaisser-Hillas function is modified into a function of one parameter only. 
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The family of one-parametric Gaisser-Hillas functions is shown in Fig. 2 (right). The prob- 
lem of finding an inverse, 

g{x; Xmax) = a (10) 

for < fl ^ 1, can be rewritten into 

- — exp ( —] = -fli/-^--- e-\ (11) 



According to the definition (1), the two (real) solutions for x are obtained from the two 
branches of the Lambert W function, 

xi,2 = -Xmax Wo,-i(-fl^/"— e-^) = -wWo,-i(- (12) 

Note that the branch — 1 or simply chooses the right or left side relative to the maximum, 
respectively. 

2 Numerics 

Before moving to the actual implementation let us review some of the possible nimerical 
and analytical approaches. 

2.1 Recursion 

For X > and W{x) > we can take the natural logarithm of (3) and rearrange it, 

W(x) = Inx -lnW(;c). (13) 

It is clear, that a possible analytical expression for W(x) exhibits a degree of self similarity. 
The W(x) function has multiple branches in the complex domain. Due to the x > and 
W(x) > conditions, the Eq. (13) represents the positive part of the Wo(x) principal branch, 
but as it turns out, in this form it is suitable for evaluation when Wo(x) > 1, i.e. when x > e. 

Unrolling the self-similarity (13) as a recursive relation, one obtains the following curi- 
ous expression for Wo(x), 

Wo(x) = Inx - \n{lnx- \n{\nx - ...)), (14) 

or in a shorthand of a continued logarithm, 

Wo(x)=ln-V- (15) 

The above expression is clearly a form of successive approximation, the final result given 
by the limit, when it exists. 

For X < and W(x) < we can multiply both sides of Eq. (3) with —1, take logarithm, 
and rewrite it to get a similar expansion for the W_i(x) branch, 

W(x) = In(-x) -ln(-W(x)). (16) 

Again, this leads to a similar recursive expression. 



W_i(x) = In(-x) - ln(-(ln(-x) - ln(-(ln(-x) - . . . )))), (17) 



or as a continued logarithm. 



W-i(x)=ln ^ . (18) 



For this continued logarithm we will use the symbol R^|(x) where n denotes the depth of 
the recursion. 

Starting from yet another rearrangement of Eq. (3), 

W(^) = SttT' (19) 

expW(x) 

we can obtain a recursion relation for the —e^^ < x < e part of the principal branch 

Wo(x) < 1, 

W.(.) ^ ,20, 

2.2 Halley's iteration 

We can apply Halley's root-finding method [8] to derive an iteration scheme for f{y) = 
W(i/) — X by writing the second-order Taylor series 

fiy) = fivn) +f'{yn) {y - Vn) + \f\yn) {y - ynf + • ■ ■ (21) 

Since root y of /(y) satisfies f{y) = we can approximate the left-hand side of Eq. (21) 
with and replace y with yn+i- Rewriting the obtained result into 

=y„ fi}hl (22) 

f'iyn) + ^f"{yn) {y„+i-y„) 

and using Newton's method yn+i — yn = ~f{yn)/f"{yn) on the last bracket, we arrive at 
the expression for the Halley's iteration for Lambert function 

W„+i = W„+ \ , (23) 

lyi Sji Ufi 

where 

tn = W„expW„-x, (24) 

_ W„+2 
'""2(W„ + 1)' 

w„ = (W„ + l)expW„. (26) 

This method is of the third order, i.e. having W„ = W(x) + 0{e) will give W^+i = 
W(x) + O(e^). Supplying this iteration with sufficiently accurate first approximation of the 
order of 0(10^^) will thus give a machine-size floating point precission O{10^^^) in at least 
two iterations. 



2.3 Fritsch's iteration 



For both branches of Lambert function a more efficient iteration scheme exists [9], 

W„+i = W„(l + £„), (27) 
where is the relative difference of successive approximations at iteration n, 

en - . (28) 



The relative difference can be expressed as 



1 + W„y \qn-2Zn 



where 



z„ = In - Wn, (30) 

(j„ =2(1 + W„) (l + W„ + ?z„). (31) 

The error term in this iteration is of a fourth order, i.e. with W„ = W{x) + 0(e„) we get 
Wn+i = W(x) + 0{el). 

Supplying this iteration with a sufficiently reasonable first guess, accurate to the order of 
0(10^^), will therefore deliver machine-size floating point precission 0(10^^^) in only one 
iteration and excessive 0(10^^"^) in two! We have to find reliable first order approximation 
that can be fed into the Fritsch iteration. Due to the lively landscape of the Lambert function 
properties, the approximations will have to be found in all the particular ranges of the 
function behavior. 



3 Initial approximations 

The following section deals with finding the appropriate initial approximations in the whole 
definition ranges of the two branches of the Lambert function. 

3.1 Branch-point expansion 

The inverse of the Lambert function, W^^ (y) = y exp y, has two extrema located at W^^ (— 1) = 
—e^^ and W^^(— oo) = 0^. Expanding W^^(i/) to the second order around the minimum 
at y = — 1 we obtain 

W-^(y)--l + ^^. (32) 

The inverse W^^ (y) is thus in the lowest order approximated with a parabolic term imply- 
ing that the Lambert function will have square-root behavior in the vicinity of the branch 
point X = —e^^, 

W_i,o(x) ~ -lTy2(l + ex). (33) 

To obtain the additional terms in expression (33) we proceed by defining an inverse func- 
tion, centered and rescaled around the minimum, i.e. f{y) = 2(eW^^(i/ — 1) + 1). Due 
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Figure 3: Successive orders of the branch-point expansion for the W_i(x) on the left 
and Wo(:t) on the right. 

to the centering and rescaling the Taylor series of this function around y = becomes 
particularly neat, 

f{y) = y^ + ly^ + \y^ + hy' + --- (34) 

Using this Taylor expansion we can derive coefficients [10] of the corresponding inverse 
function 

ri(z) = l+w(^) = (35) 
= zi/^-lz+llz3/^-i,z^ + ... (36) 

From Eq. (35) we see that z = 2(1 + Using p±{x) = ±Y^2(l + ex) as independent 
variable we can write this series expansion as 

W_i,o(x) ^ B%{x) = f^h,v^{x), (37) 

!=0 

where the lowest few coefficients hi are 



i 





1 


2 


3 


4 


5 


6 


7 


8 


9 


h 


-1 


1 


1 


11 


43 


769 


221 


680863 


1963 


226287557 


3 


72 


540 


17280 


8505 


43545600 


204 120 


37623398400 



3.2 Asymptotic series 

Another useful tool is the asymptotic expansion [2] where using 

A(fl,&) =fl-& + X]X]Qm«"'^-"'"ir+i (38) 

k tn 

where C)t^ are related to the Stirling number of the first kind, the Lambert function can be 
expressed as 

W_i,o(x) = A(ln(Tx),ln(Tln(Tx))) (39) 



with a = \nx, b = \n\nx for the Wq branch and a = ln(— x), b = ln(— ln(— x)) for the W_i 
branch. The first few terms are 

, b b(-2 + b) b(6-9b + 2b^) 

b{-12 + 36b-22b^ + 3b^) & (60 - 300b + 350b^ - I25b^ + 12b^) 
^ 12^4 + ^5 +■■■ 

3.3 Rational fits 

A useful quick-and-dirty approach to the functional approximation is to generate large 
enough sample of data points {wjexpiVi, Wj}. These points evidently lie on the Lambert 
function. Within some appropriately chosen range of values the points are fitted with a 
rational approximation 

= Hfi' 

varying the order of the polynomials in the nominator and denominator, and choosing the 
one that has the lowest maximal absolute residual in a particular interval of interest. 

For the Wo(x) branch, the first set of equally-spaced Wi component was chosen in a 
range that produced Wj exp Wi values in an interval [—0.3, 0]. The optimal rational fit turned 
out to be 

1 + flix + a2X^ + fl3X^ + a^x^ 

^O(X) -X^^^^^^ ^2^2 + b3x3 + &4X4 

where the coefficients for this first approximation Qq^' (x) are 

i 12 3 4 

Ui 5.931375839364438 11.39220550532913 7.33888339911111 0.653449016991959 
bi 6.931373689597704 16.82349461388016 16.43072324143226 5.115235195211697 



For the second fit of the Wq (x) branch a Wj range was chosen so that the Wi exp Wi values 

were produced in the interval [0.3, 2e] giving rise to the second optimal rational fit Qg ' (x) 
of the same form as in Eq. (42) but with coefficients 



Ui 2.445053070726557 1.343664225958226 0.148440055397592 0.0008047501729130 
bi 3.444708986486002 3.292489857371952 0.916460018803122 0.0530686404483322 

For the W_i(x) branch one rational approximation of the form 

^ / ^ gp + flix + aix^ 

l + hix + b2X^ + b3X^ + b^x^ + b5x5 ^ ' 

with the coefficients 
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Figure 4: Combining different approximations of Wo(x) into final piecewise func- 
tion. The number of accurate decimal places A(;c) is shown for two ranges, linear 
interval [—e^^, 0.3] on the left and logarithmic interval [0.3, 10^] on the right. The 
approximation are branch-point expansion B^^ (x) from Eq. (37) in blue, rational fits 

Qq^' (x) and Qq^' {x) from Eq. (42) in black and red, respectively, and asymptotic series 
Ao{x) from Eq. (40) in green. 



is enough. 



4 Implementation 

To quantify the accuracy of a particular approximation W{x) of the Lambert function W(x) 
we can introduce a quantity A(x) defined as 

A(x) = -log,o|W(x)-W(x)|, (44) 

so that it gives us a number of correct decimal places the approximation W(x) is producing 
for some parameter x. 

In Fig. 4 all mentioned approximations for the Wo(x) are shown in the linear interval 
[—e^^, 0.3] on the left and logarithmic interval [0.3, 10^] on the right. For each of the ap- 
proximations an use interval is chosen so that the number of accurate decimal places is 
maximized over the whole definition range. For the Wo(x) branch the resulting piecewise 
approximation 



Wo(x) 



B^Q^ (x) ; -g-i ^ X < -0.32358170806015724 

Q^Q^ (x) ; -0.32358170806015724 ^ x < 0.14546954290661823 

Q|f^ (x) ; 0.14546954290661823 ^ x < 8.706658967856612 

Ao(x) ; 8.706658967856612 ^ X < 00 



is accurate in the definition range [—e ^, 7] to at least 5 decimal places and to at least 3 
decimal places in the whole definition range. The B^^ (x) is from Eq. (37), q\^^ (x) and 
Q|f^ (x) are from Eq. (42), and Ao(x) is from Eq. (40). 



20 



15 



< 10 



p.„.,,^- ..1" ^ 









20 



15 



<1 10 




-0.3 -0.2 -0.1 0.0 0.1 0.2 0.3 

X 



1 10 100 1000 10* io5 

X 



Figure 5: Final values of the combined approximation Wo(x) (black) from Fig. 4 after 
one step of the Halley iteration (red) and one step of the Fritsch iteration (blue). 
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Figure 6: Left: Approximations of the W_i(x) branch. The branch point expansion 
Bf^[(x) is shown in blue, the rational approximation Q-i(x) in black, and the loga- 

rithmic recursion R_-^ in red. Right: Combined approximation is accurate to at least 
5 decimal places in the whole definition range. The results after applying one step 
of the Halley iteration are shown in red and after one step of the Fritsch iteration in 
blue. 



The final piecewise approximation Wo(x) is shown in Fig. 5 in black line. Using this 
approximation a single step of the Halley iteration (in red) and the Fritsch iteration (in 
blue) is performed and the resulting number of accurate decimal places is shown. As can 
be seen both iterations produce machine-size accurate floating point numbers in the whole 
definition interval except for the [9, 110] interval where the Halley method requires another 
step of the iteration. For that reason we have decided to use only (one step of) the Fritsch 
iteration in the C++ implementation of the Lambert function. 

In Fig. 6 (left) the same procedure is shown for the W_i(x) branch. The final approxi- 
mation 



W_i(x) 



'B^^\{x) ; -g-i ^ X < -0.30298541769 

Q-i{x) ; -0.30298541769 ^ x < -0.051012917658221676 (46) 
^r[^((x) ; -0.051012917658221676 ^ x < 



is accurate to at least 5 decimal places in the whole definition range [—e ^, 0] and where 
b'^J(x) is from Eq. (37), Q_i(x) is from Eq. (43), and R^^\{x) is from Eq. (18). 

In Fig. 6 (right) the combined approximation W_i (x) is shown in black line. The values 
after one step of the Halley iteration are shown in red and after one step of the Fritsch 
iteration in blue. Similarly as for the previous branch, the Fritsch iteration is superior, 
yielding machine-size accurate results in the whole definition range, while the Halley is 
accurate to at least 13 decimal places. 
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A Implementation in C++ 

Sources are available also from http : //www . ung . si/*darko/LambertW . tar . gz 



A.l Lamibert .h 

/* 

Implementation of Lambert W function 

Copyright (C) 2009 Darko Veberic, darko . vebericSung . si 

This program is free software: you can redistribute it and/or modify 
it under the terms of the GNU General Public License as published by 
the Free Software Foundation, either version 3 of the License, or 
(at your option) any later version. 

This progrEim is distributed in the hope that it will be useful, 
but WITHOUT ANY WARRANTY; without even the implied warranty of 
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 
GNU General Public License for more details. 

You should have received a copy of the GNU General Public License 
along with this program. If not, see <http://www.gnu.org/licenses/>. 

25 Jun 2009 

*/ 

#ifndef .utl .Lambert W.h_ 
#define _utl_LambertW_h_ 



/** Approximate Lambert W function 

Accuracy at least 5 decimal places in all definition range. 
See LambertWC) for details. 

\param branch: valid values are and -1 
\param x: real-valued argument \f $\geq-l/e\f $ 
\ingroup math 

*/ 

template<int branch> 

double LambertWApproximationCconst double x) ; 

/** Lambert W function 

Lambert function \f$y=-[\rm W}Cx)\f$ is defined as a solution 
to the \f$x=ye"y\f$ expression and is also known as 
"product logarithm". Since the inverse of \f$ye"y\f$ is not 
single-valued, the Lambert function has two real branches 
\f${\rm W>_0\f$ and \f${\rm W}^-[-l>\f$. 

\f${\rm W}_OCx)\f$ has real values in the interval 

\f$[-l/e,\infty]\f$ and \f$-[\rm W}_-[-l>Cx) \f $ has real values 

in the interval \f $ [-1/e , 0] \f $ . 

Accuracy is the nominal double type resolution 

(16 decimal places) , 

\param branch: valid values are and -1 

\param x: real-valued argument \f $\geq-l/e\f $ (range depends on 
branch) 

*/ 

template<int branch> 

double LambertW(const double x) ; 



#endif 



A.2 LcLmbert . cc 

/* 

Implementation of Lcimbert W function 

Copyright (C) 2009 Darko Veberic, darko . vebericSung . si 

This program is free software: you can redistribute it and/or modify 
it under the terms of the GNU General Public License as published by 
the Free Software Foundation, either version 3 of the License, or 
(at your option) any later version. 

This program is distributed in the hope that it will be useful, 
but WITHOUT ANY WARRANTY; without even the implied warranty of 
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 



GNU General Public License for more details. 

You should have received a copy of the GNU General Public License 
along with this program. If not, see <http://www.gnu.org/licenses/>. 

*/ 

#include <iostream> 
#include <cmath> 
#include <limits> 
#include "LambertW. h" 

using namespace std; 
namespace Lamb ertWDe tail { 
const double kInvE = 1/M_E; 



template<int n> 

inline double BranchPointPolynomiaKconst double p) ; 



templateO 

inline 

double 

BranchPointPolynomial<l>Cconst double p) 
{ 

return 

-1 + p; 

} 



templateO 

inline 

double 

BranchPointPolynomial<2> (const double p) 
{ 

return 

-1 + p*(l + p*(-l./3)); 

} 



templateO 

inline 

double 

BranchPointPolynomial<3> (const double p) 
{ 

return 

-1 + p*(l + p*(-l./3 + p*(ll./72))); 

> 



templateO 

inline 

double 

BranchPointPolynomial<4> (const double p) 
{ 

return 

-1 + p*(l + p*(-l./3 + p*(ll./72 + p*(-43./540)))) ; 

} 



templateO 

inline 

double 

BranchPointPolynomial<5> (const double p) 
{ 

return 

-1 + p*(l + p*(-l./3 + p*(ll./72 + p*(-43./540 + 
p*(769./17280))))); 

} 



templateO 

inline 

double 

BranchPointPolynomial<6> (const double p) 
{ 

return 

-1 + p*(l + p*(-l./3 + p*(ll./72 + p+(-43./540 + p* (769 . /17280 
+ P+C-221./8505)))))) ; 



> 



templateO 

Inline 

double 

BranchPointPolynomial<7> (const double p) 

{ 

return 

-1 + p«(l + p*(-l./3 + p*(ll./72 + p»(-43./540 + p*(769./17280 
+ p*(-221./8606 + p*(680863./43646600)) ))))); 

> 



templateO 

inline 

double 

BranchPointPolynomial<8> (const double p) 

{ 

return 

-1 + p*(l + p*C-l./3 + p*(ll./72 + p*(-43./540 + p*(769./17280 
+ p* (-221 . /8505 + p* (680863. /43545600 + p*(-1963./204120) ))))))) ; 

} 



templateO 

inline 

double 

BranchPointPolynomial<9> (const double p) 

{ 

return 

-1 + p»(l + p*(-l./3 + p»(ll./72 + p»C-43./B40 + p* (769 . /17280 
+ p*(-221./8505 + p* (680863. /43545600 + p* (-1963 . /204120 
+ p» (226287667 . /37623398400 .))))))))); 

> 



template<int order> 

inline double AsymptoticExpaQsion(const double a, const double b) ; 



templateO 

inline 

double 

AsymptoticExpansion<0>(const double a, const double b) 

{ 

return a - b; 

> 



templateO 

inline 

double 

AsymptoticExpansion<l>(const double a, const double b) 

{ 

return a - b + b / a; 

> 



templateO 

inline 

double 

AsymptoticExpansion<2>(const double a, const double b) 

{ 

const double ia = 1 / a; 

return a-b+b/a* (1+ia* 0.5*(-2 + b)); 

} 



templateO 

inline 

double 

AsymptoticExpansion<3>(const double a, const double b) 

{ 

const double ia = 1 / a; 
return a-b+b/a* 
(1+ia * 

(0.6*(-2 + b) + ia * 

1/6.* (6 + b*(-9 + b*2)) 

) 

); 

> 



templateO 

inline 

double 

AsyinptoticExpansion<4>(const double a, const double b) 
{ 

const double ia = 1 / a; 



return a-b+b/ a* 
(1+ia * 

(0.6*(-2 + b) + ia * 

(1/6. *(6 + b»(-9 + b*2)) + ia * 

1/12.*(-12 + b*(36 + b*(-22 + b*3))) 

) 

) 

)i 

} 



templateO 

inline 

double 

AsymptoticExpansion<5> (const double a, const double b) 

< 

const double ia = 1 / a; 
return a-b + b/ a* 
(1 + ia » 

(0.5*(-2 + b) + ia » 

(1/6. »(6 + b*(-9 + b»2)) + ia * 

(1/12.*(-12 + b»(36 + b*(-22 + b*3))) + ia » 

1/60. *(60 + b*(-300 + b*(350 + b*(-125 + b*12)))) 

) 

) 

) 

)i 

} 

template<int branch> 
class Branch -C 

public: 

template<int order> 

static double BranchPointExpansion(const double x) 

{ return BranchPointPolynomial<order>(eSign * sqrt(2*CM_E*x+l))) ; } 

// Asymptotic expansion 

// Corless et al. 1996, de Bruijn (1981) 

template<int order> 

static 

double 

AsymptoticExpansion(const double x) 
{ 

const double logsx = log(eSign * x) ; 

const double logslogsx = log(eSign * logsx); 

return Lambert UDetail : :AsymptoticExpansion<order>(logsx, logslogsx) ; 

} 

tGmplate<int n> 

static inline double RationalApproximation(const double x) ; 

// Logarithmic recursion 
template<int order> 

static inline double LogRecursion(const double x) 
{ return LogRecursionStep<order>(log(eSign * x) ) ; } 

// generic approximation valid to at least 5 decimal places 
static inline double Approximation(const double x) ; 

private : 

// enum { eSign = 2*branch + 1 }; this doesn't work on gcc 3.3.3 
static const int eSign = 2*branch + 1; 

template<int n> 

static inline double LogRecursionStep(const double logsx) 

■[ return logsx - log(eSign * LogRecursionStep<n-l> (logsx) ) ; } 

>; 



// Rational approximations 

templateO 

templateO 

inline 

double 

Branch<0>: :RationalApproximation<0>(const double x) 

{ 

return 3:*(60 + x*(114 + 17*i)) / (60 + 3t*(174 + 101*3t))i 

} 



templateO 
templateO 
inline 
double 

Branch<0>; :RationalApprcximation<l>(const double x) 
{ 

// branch 0, valid for [-0.31,0.3] 



return 

X * (1 + X * 

(6.931375839364438 + x * 
(11.392205505329132 + x « 

(7.338883399111118 + x*0. 6634490169919699) 

) 

) 

) / 

(1 + X * 

(6.931373689597704 + x • 
(16.82349461388016 + x » 

(16.43072324143226 + x*6. 116236196211697) 

) 

) 

); 

> 



templateO 
templateO 
inline 
double 

Branch<0>: :RationalApproximation<2>(const double x) 

{ 

// branch 0, valid for [-0.31,0.6] 
return 

X * (1 + X * 

(4.790423028527326 + x * 

(6.696946076293267 + x * 2.4243096806908033) 

) 

) / 

(1 + X * 

(6.790432723810737 + x • 
(10.986446930034288 + x * 

(7.391303898769326 + x • 1.1414723648617864) 

) 

) 

); 

> 



templateO 
templateO 
inline 
double 

Branch<0>: :RationalApproximation<3>(const double x) 

{ 

// branch 0, valid for [0.3,7] 
return 

X * (1 + X * 

(2.4460530707265568 + x * 
(1.3436642269682266 + x * 

(0.14844006639769196 + x * 0.0008047601729129999) 

) 

) 

) / 

(1 + X » 

(3.4447089864860025 + X * 
(3.2924898573719623 + x * 

(0.9164600188031222 + x * 0.06306864044833221) 

) 

) 

); 

} 



templateO 
templateO 
inline 
double 

Braiich<-1>: :RationalApproximation<4> (const double x) 

{ 

// branch -1, valid for [-0.3,-0.06] 
return 

(-7.814176723907436 + x » 

(263.88810188892484 + x » 667.9493176902304) 
) / 

(1 + I • 

(-60.43968713690808 + x * 
(99.98567083107612 + x * 
(682.6073999909428 + x » 

(962.1784396969866 + i * 1477.9341280760887) 

) 

) 

) 

); 

} 



// stopping conditions 



templateO 
templateO 
inline 
double 

Branch<0>: :LogRecursionStep<0> (const double logsx) 

{ 

return logsx; 

} 



templateO 
templateO 
inline 
double 

Branch<-1>: :LogRecursionStep<0>(const double logsx) 

{ 

return logsx; 

} 



templateO 

inline 

double 

Branch<0> : : Approximation(const double x) 
{ 

if Cx < -0.32358170806015724) { 
if (x < -kInvE) 

return numeric_limits<double> : :quiet_NaN() ; 
else if (x < -kInvE+le-5) 

return BranchPointExpaiision<5> (x) ; 
else 

return BranchPointExpansion<9> (x) ; 
} else { 

if (x < 0.14546954290661823) 

return RationalApproximation<l>(x) ; 
else if (x < 8.706658967856612) 

return RationalApproximation<3>(x) ; 
else 

return AsymptoticExpansion<5>(x) ; 

> 



templateO 

inline 

double 

Branch<-1>: :Approximation(const double x) 
{ 

if Cx < -0.051012917658221676) { 
if (x < -kInvE+le-5) { 
if (x < -kInvE) 

return numeric_limits<double> : :quiet_NaN() ; 
else 

return BranchPointExpansion<5>(x) ; 
} else { 

if (x < -0.30298541769) 

return BranchPointExpansion<9>(x) ; 
else 

return RationalApproximation<4>(x) ; 

> 

} else { 
if (x < 0) 

return LogRecursion<9> (x) ; 
else if (x == 0) 

return -numeric_limits<double> : :infinity{) ; 
else 

return numeric_limits<double> : :quiet_MaN() ; 

> 



// iterations 

inline 
double 

HalleyStepCconst double x, const double w) 

{ 

const double ew = expCw); 
const double wew = w * ew; 
const double wewx = wew - x; 
const double wl = w + 1 ; 

return w - wewx / (ew * wl - (w + 2) * wewx/(2*wl)) ; 

} 



inline 
double 

FritschStep(const double x, const double w) 

{ 



const double z = log(x/w) - w; 
const double wl = w + 1; 

const double q = 2 * wl * Cwl + (2/3. )*z); 
const double eps = z / wl * Cq - z) / (q - 2*z) ; 
return w * Cl + eps); 



tempi at e< 

double IterationStepCconst double x, const double w) 

> 

Inline 
double 

IterateCconst double x, double w, const double eps = le-6) 
{ 

for Cint i = 0; i < 100; { 

const double ww = IterationStepCx, w) ; 
if CfabsCww - w) <= eps) 

return ww; 
w = ww; 

} 

cerr « "convergence not reached." « endl; 
return w; 



template< 

double IterationStepCconst double x, const double w) 

> 

struct Iterator { 

static 
double 

DoCconst int n, const double x, const double w) 
{ 

for (int i = 0; i < n; ++i) 

w = IterationStepCx, w) ; 
return w; 

} 

tempi at e< int n> 

static 

double 

DoCconst double x, const double w) 
{ 

for Cint i = 0; i < n; ++i) 

w = IterationStepCx, w) ; 
return w; 

> 

template<int n, class = void> 

struct Depth { 

static double Recurse Cconst double x, double w) 

■[ return Depth<n-1>: iRecurseCx, IterationStepCx, w)); } 

}; 

// stop condition 
template<class T> 
struct Depth<l, T> { 

static double RecurseCconst double x, const double w) 

{ return IterationStepCx, w) ; } 



}; 

// identity 
teDiplate<class T> 
struct Depth<0, T> { 

static double RecurseCconst double x, const double w) 

{ return w; } 

}; 

}; 



} // LambertWDetail 



teDiplate<int branch> 
double 

LambertWApproximationCconst double x) 
i 

return LajnbertWDetail : : Branch<branch> : rApproximationCx) ; 

} 

// instantiations 

template double LambertWApproximation<0> Cconst double x) ; 
template double LambertWApproximation<-l> (const double x) ; 



template<int branch> 

double LambertU Cconst double x) ; 

templateo 
double 

LambertW<0> (const double x) 
{ 

if CfabsCx) > le-6 M x > -LambertWDetail: :kInvE + le-5) 
return 

LambertWDetail: : 

Iterator<LambertWDetail : :FritschStep>: : 
Depth<l>: : 

Recurse Cx , LambertWApproximation<0> Cx) ) ; 

else 

return LambertWApproximation<0>Cx) ; 

} 

tempi at e<> 
double 

LambertW<-l> Cconst double x) 
{ 

if Cx > -LambertWDetail: :kInvE + le-5) 

return 

LajnbertWDetail: : 

Iterator<LambertWDetail : :FritschStep>: : 
Depth<l>: : 

Recurse Cx , Lambert WApproximation<-l> Cx) ) ; 

else 

return LambertWApproximation<-l>Cx) ; 

} 

// instantiations 

template double LambertW<0> Cconst double x) ; 
template double LambertW<-l> Cconst double x) ; 



